family = binomial(),
data = y.fall)
comp2 <- proc.time() - start.time
newdat <- with(y.fall, expand.grid(Time = seq(0, 24, length.out = 24)))
predict_unimodal.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
predict1=ggeffects::predict_response(model=unimodal,margin="empirical",
terms="Time[all]",type = "response")
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
ggeffects::predict_response(unimodal2,terms="Time[all]",newdata=newdat)
?predict_response
unimodal2
unimodal2$frame
unimodal2 <- glmmTMB(cbind(success, failure) ~ Time,
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]",newdata=newdat)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(Time),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
y.fall2$Time
y.fall2=y.fall
y.fall2$Time
y.fall2$Cos.Time=cos(2*pi*y.fall2$Time/24)
y.fall2$Sin.Time=sin(2*pi*y.fall2$Time/24)
unimodal2 <- glmmTMB(cbind(success, failure) ~Cos.Time+Sin.Time,
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~Cos.Time+Sin.Time,
family = binomial(),
data = y.fall2)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time) + sin(2*pi*Time),
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time) + sin(2*pi*Time),
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~ (2*pi*Time/24) + (2*pi*Time/24),
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~ 2*pi*Time/24 + 2*pi*Time/24,
family = binomial(),
data = y.fall)
y.fall
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24),
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(Time/24) + sin(Time/24),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
pi
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24)+
(1 + cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24) || Site),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24)+
(1 + cos(2*pi*Time/24) + sin(2*pi*Time/24) || Site),
family = binomial(),
data = y.fall)
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24)+
(1 + cos(2*pi*Time/24) + sin(2*pi*Time/24) || Site),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
ggeffects::predict_response(unimodal2,terms="Time[all]",condition = (pi=3.141593))
?pi
unimodal <- mixed_model(fixed = cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24),
random = ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24) || Site,
family = binomial(),
data = y.fall)
predict1=ggeffects::predict_response(model=unimodal,margin="empirical",
terms="Time[all]",type = "re")
predict1=ggeffects::predict_response(model=unimodal,terms="Time[all]")
unimodal
predict1=ggeffects::predict_response(model=unimodal)
predict1=ggeffects::predict_response(model=unimodal,terms="Time[all]")
predict1=ggeffects::predict_response(model=unimodal,margin="empirical",
terms="Time[all]")
predict1=ggeffects::predict_response(model=unimodal,margin="empirical",
terms="Time[all]",type = "random")
predict1=ggeffects::predict_response(model=unimodal,margin="empirical",
terms="Time[all]",type = ".")
predict1=ggeffects::predict_response(model=unimodal,margin="mean_reference",
terms="Time[all]",type = "random")
predict1=ggeffects::predict_response(model=unimodal,margin="mean_reference",
terms="Time[all]")
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24)+
(1 + cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24) || Site),
family = binomial(),
data = y.fall)
ggeffects::predict_response(unimodal2,terms="Time[all]")
unimodal <- mixed_model(fixed = cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24),
random = ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24) || Site,
family = binomial(),
data = y.fall)
predict1=ggeffects::predict_response(model=unimodal,terms="Time[all]")
ggeffects::predict_response(unimodal2,terms="Time[all]"
predict1=ggeffects::predict_response(model=unimodal,terms="Time[all]")
ggeffects::predict_response(unimodal2,terms="Time[all]")
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]")
predicted2
print(predicted2)
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2)
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]")
print(predicted2,n=24)
y.fall
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", condition = c(Site="10A"))
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", condition = c(Site="10A"),type="random")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", type="random")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", type="random",back_transform=FALSE)
print(predicted2,n=24)
unimodal2
summary(unimodal2)
predict(unimodal2)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", type="random",back_transform=FALSE,margin="mean_reference")
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", type="random",back_transform=FALSE,margin="mean_reference")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]", type="random",back_transform=FALSE,margin="marginalmeans")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]",back_transform=FALSE,margin="marginalmeans")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]",back_transform=FALSE,margin="empirical")
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]",back_transform=FALSE,margin="marginalmeans")
predicted2
print(predicted2,n=24)
predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]",back_transform=TRUE,margin="marginalmeans")
print(predicted2,n=24)
Time = seq(0, 24, length.out = 24)
Time = seq(0, 24, length.out = 24)
newdat2=data.frame(Time = seq(0, 24, length.out = 24),Site=NA)
predict(unimodal2,newdata = newdat2,type="response")
predict(unimodal2,newdata = newdat,type="response")
newdat <- with(y.fall, expand.grid(Time = seq(0, 24, length.out = 24),Site=NA))
predict_unimodal.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
predict_unimodal.marginal
predict_unimodal.marginal[,3:5]
predict_unimodal.marginal[,3:5]=plogis(predict_unimodal.marginal[,3:5])
predict_unimodal.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*pi*Time/24) + sin(2*pi*Time/24)+
(1 + cos(2*pi*Time/24) + sin(2*pi*Time/24) || Site),
family = binomial(),
data = y.fall)
predict2_unimodal.marginal=predict(unimodal2,newdata = newdat,type="response")
predict2_unimodal.marginal
predict_unimodal.marginal[,2]
predict_unimodal.marginal
plogis(predict_unimodal.marginal[,3])
predict_adaptive.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
predict.TMB.marginal=predict(unimodal2,newdata = newdat,type="link",se.fit = TRUE)
predict.TMB.marginal
predict.TMB.marginal=predict(unimodal2,newdata = newdat,type="link",se.fit = TRUE)
predict.TMB.marginal$LCL = predict.TMB.marginal$fit-1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal$UCL = predict.TMB.marginal$fit+1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal
data.frame(predict.TMB.marginal)
predict.TMB.marginal=data.frame(predict.TMB.marginal)
predict.TMB.marginal
predict.TMB.marginal=data.frame(plogis(predict.TMB.marginal[,-2]))
predict.TMB.marginal=data.frame(predict.TMB.marginal[,-2])
predict.TMB.marginal
apply(predict.TMB.marginal,2,plogis)
predict.TMB.marginal.transformed=apply(predict.TMB.marginal,2,plogis)
predict_adaptive.marginal
predict_adaptive.marginal.transormed=apply(predict_adaptive.marginal,2,plogis)
predict_adaptive.marginal.transormed
predict_adaptive.marginal.transormed=apply(predict_adaptive.marginal[,2:4],2,plogis)
predict_adaptive.marginal.transormed
predict_adaptive.marginal.transormed=apply(predict_adaptive.marginal[,3:5],2,plogis)
predict_adaptive.marginal.transormed
?effectPlotData
head(predict.TMB.marginal.transformed)
predict_adaptive.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
predict_adaptive.marginal
newdat
predict.adaptive.marginal
predict.adaptive.marginal <- effectPlotData(unimodal,newdata=newdat, marginals=TRUE)
predict.adaptive.marginal
predict.adaptive.marginal[,3:5] <- apply(predict_adaptive.marginal[,3:5],2,plogis)
predict.adaptive.marginal
predict.TMB.marginal=predict(unimodal2,newdata = newdat,type="link",se.fit = TRUE)
predict.TMB.marginal$LCL = predict.TMB.marginal$fit-1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal$UCL = predict.TMB.marginal$fit+1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal
predict.TMB.marginal=data.frame(apply(predict.TMB.marginal[,-2],2,plogis))
predict.TMB.marginal
predict.TMB.marginal=data.frame(predict.TMB.marginal)
predict.TMB.marginal=apply(predict.TMB.marginal[,-2],2,plogis)
predict.TMB.marginal
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l")
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l")
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=2,lty=3)
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.0015))
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=3)
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002))
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=3)
predict.TMB.marginal
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002))
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$fit,lwd=2,col=2,lty=3)
predict.TMB.marginal=data.frame(apply(predict.TMB.marginal[,-2],2,plogis))
predict.TMB.marginal
predict.TMB.marginal=predict(unimodal2,newdata = newdat,type="link",se.fit = TRUE)
predict.TMB.marginal$LCL = predict.TMB.marginal$fit-1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal$UCL = predict.TMB.marginal$fit+1.96*predict.TMB.marginal$se.fit
predict.TMB.marginal=data.frame(predict.TMB.marginal)
predict.TMB.marginal
predict.TMB.marginal=data.frame(apply(predict.TMB.marginal[,-2],2,plogis))
predict.TMB.marginal
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002))
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$fit,lwd=2,col=1,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$LCL,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$UCL,lwd=2,col=2,lty=3)
#unimodal2 <- glmmTMB(cbind(success, failure) ~ cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24)+
#                         (1 + cos(2*3.141593*Time/24) + sin(2*3.141593*Time/24) || Site),
#                        family = binomial(),
#                        data = y.fall)
#predicted2=ggeffects::predict_response(unimodal2,terms="Time[all]",back_transform=TRUE,margin="marginalmeans")
#print(predicted2,n=24)
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002),lty=1)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=1,lty=1)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=1)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$fit,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$LCL,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$UCL,lwd=2,col=2,lty=3)
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002),lty=1,
ylab="Predicted Probability of Activity",xlab="Time (hour)")
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=2,col=1,lty=1)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=2,col=2,lty=1)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$fit,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$LCL,lwd=2,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$UCL,lwd=2,col=2,lty=3)
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=2,ylim=c(0,0.002),lty=1,lwd=3,
ylab="Predicted Probability of Activity",xlab="Time (hour)")
plot(predict.adaptive.marginal$Time,predict.adaptive.marginal$pred,type="l",lwd=3,ylim=c(0,0.002),lty=1,
ylab="Predicted Probability of Activity",xlab="Time (hour)")
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$low,lwd=3,col=1,lty=1)
lines(predict.adaptive.marginal$Time,predict.adaptive.marginal$upp,lwd=3,col=2,lty=1)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$fit,lwd=3,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$LCL,lwd=3,col=2,lty=3)
lines(predict.adaptive.marginal$Time,predict.TMB.marginal$UCL,lwd=3,col=2,lty=3)
newdat
y.fall$Site
unique(y.fall$Site)
?with
newdat <- expand.grid(Time = seq(0, 24, length.out = 24),Site=unique(y.fall$Site))
newdat
detach("package:ggeffects", unload = TRUE)
install.packages("marginaleffects")
diurnal= c(0.2,0.4,0.8,0.8,0.8)
diurnal= c(0.2,0.4,0.8,0.8,0.8)
n=c(2,2,10,10,10)
sum(diurnal)
sum(diurnal/n)
sum(diurnal/n)/length(diurnal)
mean(diurnal)
sum(diurnal)/length(diurnal)/n
sum(diurnal)/n/length(diurnal)
sum(diurnal)^(1/n)
sum(diurnal^(1/n))
sum(diurnal^(1/n))/length(diurnal)
mean(diurnal)
sum(diurnal^(1/n))/length(diurnal)
diurnal= c(0.2,0.2,0.8,0.8,0.8)
n=c(2,2,10,10,10)
mean(diurnal)
sum(diurnal^(1/n))/length(diurnal)
mean(diurnal)
sum(diurnal^(1/n))/length(diurnal)
n=c(2,2,100,100,100)
mean(diurnal)
sum(diurnal^(1/n))/length(diurnal)
mean(diurnal)
sum(diurnal^(1/n))/length(diurnal)
n=c(2,2,100,100,100)
sum(diurnal^(1/n))/length(diurnal)
diurnal^(1/n)
weighted.mean(diurnal, w=n, na.rm = FALSE)
?weighted.mean
lm(diurnal~1, weights=n)
weighted.mean(diurnal, w=n, na.rm = FALSE)
?weighted.mean
diurnal = c(0.2,0.2,0.8,0.8,0.8)
n = c(2,2,10,10,10)
mean(diurnal)
#Good. Same
weighted.mean(diurnal, w=n, na.rm = FALSE)
lm(diurnal~1, weights=n)
install.packages("robumeta")
library(robumeta)
data(hierdat) # Small-Sample Corrections- Hierarchical Dependence Model HierModSm <- robu(formula = effectsize ~ binge + followup + sreport + age, data = hierdat, studynum = studyid, var.eff.size = var, modelweights = "HIER", small = TRUE) print(HierModSm) # Output results
var
hierdat$var
studyid
hierdat$studyid
diurnal = c(0.2,0.2,0.8,0.8,0.8)
n = c(2,2,10,10,10)
lm(diurnal~1, weights=n)
mean(diurnal) # not good
weighted.mean(diurnal, w=n, na.rm = FALSE)
lm(diurnal~1, weights=n)
sum(diurnal^(1/n))/length(diurnal)
sum(diurnal^(1/n))/length(diurnal)
?glmmTMB::glmmTMB
remotes::install_github("saudiwin/ordbetareg_pack",build_vignettes=TRUE,dependencies = TRUE)
remotes::install_github("saudiwin/ordbetareg_pack",build_vignettes=TRUE,dependencies = TRUE)
remove.packages("htmltools")
install.packages("htmltools")
install.packages("microbenchmark")
library(microbenchmark)
set.seed(1)
m <- 1e4
n <- 1e3
k <- 3e2
X <- matrix(rnorm(m*k), nrow=m)
Y <- matrix(rnorm(n*k), ncol=n)
microbenchmark({
Z <- X %*% Y
}, times=10L)
y = 3
install.packages('unmarked')
object = seq(1,10)
object
matrix = matrix(seq(1,10),nrow=2,ncol=5)
matrix[2,]
install.packages("styler")
mtcars
mtcars[,1]
Which(mtcars$mpg>20)
which(mtcars$mpg>20)
mtcars$mpg[which(mtcars$mpg>20)] = NA
mtcars
mtcars$mpg
mtcars$mpg
mean(mtcars$mpg,na.rm=TRUE)
mean(mtcars$mpg,
na.rm=TRUE
)
sd(mean(mtcars$mpg,na.rm=TRUE)
)
sd(mtcars$mpg,na.rm=TRUE)
sim = rnorm(mean = mean(mtcars$mpg,na.rm=TRUE),
sd = sd(mtcars$mpg,na.rm=TRUE)
)
# Wrapping functions
sim = rnorm(10,
mean = mean(mtcars$mpg,na.rm=TRUE),
sd = sd(mtcars$mpg,na.rm=TRUE)
)
sim
sim
mu_baseline = 12.5
mu_baseline = 12.5
sd_bseline = 1
x = rnorm(1000, mu_baseline,sd_bseline)
x
mu_baseline+1.5*sd_bseline
x+1.5*sd_bseline
hist(x+1.5*sd_bseline)
x2=x+1.5*sd_bseline
hist(x2)
mean(x2)
sd(x2)
install.packages("pwr")
library(pwr)
?pwr
setwd("~/GitHub/bgerber123.github.io/classfiles/bayeslab/Custom MCMC Example")
####  Source the function
source("binom.logit.mcmc.R")
####  Prior specification
alpha=1
beta=1
#Plot of prior
curve(dbeta(x,shape1=alpha,shape2=beta),lwd=3,
xlab="Probability",ylab="Probabilty Density",
main="Prior Probability of Success",ylim=c(0,2))
#Plot of prior
curve(dbeta(x,shape1=alpha,shape2=beta),lwd=3,
xlab="Probability",ylab="Probabilty Density",
main="Prior Probability of Success",ylim=c(0,2))
y=c(0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0)
#How many iterations to run the algorithm
n.mcmc=1000
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=0.2,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
#How many iterations to run the algorithm
n.mcmc=5000
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=0.2,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.theta/binom.logit.out$n.mcmc))
hist(binom.logit.out$p.save,breaks=30,xlim=c(0,1),prob=TRUE,col=8,
main=expression(paste("Posterior for ",p)),xlab=expression(p))
# remove the first 500 iterations as 'burn-in'
mean(binom.logit.out$p.save[-(1:500)])
quantile(binom.logit.out$p.save[-(1:500)],c(0.025,0.975))
#posterior median
quantile(binom.logit.out$p.save[-(1:500)],probs=0.5)
# 95% Creditble Intervals
quantile(binom.logit.out$p.save[-(1:500)],c(0.025,0.975))
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=0.4,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=0.7,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=0.8,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.theta/binom.logit.out$n.mcmc))
hist(binom.logit.out$p.save,breaks=30,xlim=c(0,1),prob=TRUE,col=8,
main=expression(paste("Posterior for ",p)),xlab=expression(p))
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.theta/binom.logit.out$n.mcmc))
binom.logit.out$mh.theta
binom.logit.out
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.p/binom.logit.out$n.mcmc))
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=1,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.p/binom.logit.out$n.mcmc))
hist(binom.logit.out$p.save,breaks=30,xlim=c(0,1),prob=TRUE,col=8,
main=expression(paste("Posterior for ",p)),xlab=expression(p))
#How many iterations to run the algorithm
n.mcmc=10000
binom.logit.out=binom.logit.mcmc(y,
alpha,
beta,
p.tune=1,
n.mcmc=n.mcmc)
#What proportion of guesses were accepted?
#Want this b/w 0.2 an 0.6. Just not near 0 or 1
binom.logit.out$mh.p/n.mcmc
plot(binom.logit.out$p.save,type="l",ylab=expression(p),
main=paste("Acceptance: ",binom.logit.out$mh.p/binom.logit.out$n.mcmc))
hist(binom.logit.out$p.save,breaks=30,xlim=c(0,1),prob=TRUE,col=8,
main=expression(paste("Posterior for ",p)),xlab=expression(p))
#posterior mean
mean(binom.logit.out$p.save[-(1:500)])
#posterior median
quantile(binom.logit.out$p.save[-(1:500)],probs=0.5)
# 95% Credible Intervals
quantile(binom.logit.out$p.save[-(1:500)],c(0.025,0.975))
